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We investigate the quantum many-body dynamics of dissociation of a Bose-Einstein condensate 
of molecular dimers into pairs of constituent bosonic atoms and analyze the resulting atom-atom 
correlations. The quantum fields of both the molecules and atoms are simulated from first principles 
in three dimensions using the positive-P representation method. This allows us to provide an exact 
treatment of the molecular field depletion and s-wave scattering interactions between the particles, 
as well as to extend the analysis to nonuniform systems. In the simplest uniform case, we find 
that the major source of atom-atom decorrelation is atom-atom recombination which produces 
molecules outside the initially occupied condensate mode. The unwanted molecules are formed 
from dissociated atom pairs with nonopposite momenta. The net effect of this process - which 
becomes increasingly significant for dissociation durations corresponding to more than about 40% 
conversion - is to reduce the atom-atom correlations. In addition, for nonuniform systems we find 
that mode mixing due to inhomogeneity can result in further degradation of the correlation signal. 
We characterize the correlation strength via the degree of squeezing of particle number-difference 
fluctuations in a certain momentum-space volume and show that the correlation strength can be 
increased if the signals are binned into larger counting volumes. 
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I. INTRODUCTION 

Dissociation of a diatomic molecule produces two 
quantum mechanically entangled atoms with equal and 
opposite momenta in the molecule's rest frame. These 
atoms have Einstein-Podolsky-Rosen (EPR) type corre- 
lations in position and momentum JJ, and hence are 
of fundamental interest 0, Q- Experimental advances 
in coherent manipulation of quantum gases are now en- 
abling the production of ultracold molecular gases Q 
and even molecular Bose-Einstein condensates (BECs) 
from atomic condensates and degenerate Fermi gases. 
These molecules can in turn be dissociated laB B into 
strongly correlated ensembles of bosonic |3, flOl llll | or 
fermionic [T^ [l3l| atoms, thus extending possible fun- 
damental tests of quantum mechanics into macroscopic 
regimes This has close analogies with experiments 
in quantum optics using continuous-variable quadrature 
correlations generated via parametric down-conversion 




Apart from the quantum-atom optical aspect of these 
and related studies (see^g_,^El iH IH Ha HI HI HI 
El m m m m Hi ES molecular dissociation 

can serve as a probe of two-body interactions, includ- 
ing collisional resonances and spectroscopic properties of 
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Feshbach resonance molecules [s^, 113, 13 Hll • 

An important recent development in the study of ultra- 
cold quantum gases is the direct measurement of atom- 
atom correlation functions and quantum statistics via the 
analysis of the noise in absorption images H, |3^ [s^, Issll . 
atom counting using microchannel plate detectors p39l 
Eoj , and fluorescence imaging combined with high- finesse 
optical cavities 0, 0. These experiments have close 
parallels with the pioneering photon correlation mea- 
surements of Hanbury Brown and Twiss (HBT) "43] and 
initiated further theoretical studies of HBT interferome- 
try as a sensitive tool for probin g qu antum many-body 
states of ultracold gases 0, EElM Ell- All these 
correlation measurements are ultimately related to the 
quantum-atom optics counterpart of Glauber's second- 
order, or density-density, correlation function While 
giving access to the simplest higher-order correlations, 
these techniques may in the future lead to more sophisti- 
cated correlation measurements, such as between matter- 
wave quadratures required to demonstrate the version of 
the EPR paradox proposed in Ref. 0| . 

In correlation measurements via the noise in absorp- 
tion images of ballistically expanded clouds H , the 
extracted correlation signal reflects the momentum corre- 
lations before expansion ^37j . This is particularly suitable 
for studies of molecular dissociation, since the correla- 
tions arise from momentum conservation and are there- 
fore most fundamentally between atoms with equal and 
opposite momenta, Kk and — ftk. 

In this paper, we study the quantum dynamics of dis- 
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sociation of a molecular BEC and analyze the atom- 
atom correlations in momentum space. Spatial corre- 
lations after expansion will be studied in a subsequent 
paper |5C| . Our analysis is based on first- principles 
quantum simulations of the coupled atomic-molecular 
system in three dimensions (31)) usin g th e positive P- 
representation method ISZ . Is^ . Is^. l55l| . This allows 
us to go beyond the analytical results known for the case 
of a spatially uniform condensate and undepleted, classi- 
cal molecular field 0, , as well as beyond the pairing 
mean- field method used in Ref. |l2l |. In previous re- 
lated work our simulations were limited to one spatial 
dimension and to the calculation of the number differ- 
ence squeezing between all atoms with positive momenta 
and all with negative momenta finl |. In this paper we 
simulate three spatial dimensions and calculate squeez- 
ing localized in momentum space. 

In the present treatment, we are able to extend these 
results and analyze the molecular depletion without ap- 
proximations. In addition, the positive-P method allows 
us to treat nonuniform condensates, as well as s-wave 
scattering interactions between the species. The main 
limitation of the method is that the positive-P simula- 
tions eventually fail as the simulation time increases, par- 
ticularly when two-body scattering is included |5^. For 
this reason, the simulations involving s-wave interactions 
are performed with slightly reduced values of the scatter- 
ing lengths, or else are limited to short dissociation times. 
Nevertheless, the results allow us to gain quantitative in- 
sights into the role of these interactions in atom-atom 
correlations, which have not been analyzed before. 

More importantly, we find that mode mixing due to 
inhomogeneity of the condensates can potentially have 
a quite disruptive effect on the atom- atom correlations. 
We give quantitative estimates of this effect. Similarly, 
we quantify the disruptive role of atom-atom recombi- 
nation which becomes increasingly important past the 
initial spontaneous dissociation regime. 

We emphasize that for the correct interpretation of 
the results on correlations, especially when comparing 
the uniform and nonuniform results, it is important to 
operate with observables that have a well-defined opera- 
tional meaning. We find that the correlation strength is 
most conveniently quantified via the normalized number- 
difference fluctuations. This is defined with respect to the 
number of atoms in a certain momentum-space volume 
around pairs of the carrier momenta of interest. This 
measure of the correlation strength can be defined with 
respect to different counting volumes, which corresponds 
to the procedure of binning often employed in experi- 
ments. We find that the strength of the correlation signal 
increases for larger bin sizes, which is in agreement with 
the results of analogous experimental measurements of 
Ref. 13 performed in position space. 

This paper is organized as follows. SectionlTlldescribes 
the system we study including relevant approximations. 
The positive-P simulation method and the full stochastic 
equations we simulate are presented. Section [llll reviews 



the idealized analytic model presented by Kherutsyan 
[Tsf. This is the reference against which we can deter- 
mine the effect of the various realistic processes included 
in our simulations. Section IIVI is the heart of the pa- 
per and reports numerical simulations that quantify the 
degradation of quantum correlations in the presence of: 
molecular depletion (Sec. IIV B|) . scattering (Sec. IIV Cp . 
and spatial nonuniformity (Sec. IIVD|I . The effect of bin- 
ning data is briefly discussed in Sec. IIV D 31 The paper 
ends with a summary, followed by appendices on Fourier 
transforming the fields. 

II. THE MODEL 

The quantum field theory effective Hamiltonian de- 
scribing the system we shall study is given by 0, [s^l : 

The molecular and atomic fields are, respectively, de- 
scribed by the bosonic operators *I'o(x, t) and ^'i(x, i) 
satisfying the following equal-time commutation relation: 
[*.(x,t),*j(x',t)] = 5,,5(x-x') (i,j = 0,1). The first 
term in the Hamiltonian describes the kinetic energy, 
where m\ and mo = 2mi are the atomic and molecu- 
lar masses. The trapping potentials, including internal 
energies are given by Vi(x). This determines the de- 
tuning 2A which corresponds to the overall energy mis- 
match 2ftA = W(lVx{^) - Vb(0)] = 2Ey - Eo between the 
free two-atom state at the dissociation threshold and the 
bound molecular state. 

The coupling term x is responsible for coherent con- 
version of molecules into atom pairs, e.g., via Raman 
transitions or a Feshbach resonance (see, for example, 
Refs. 113 m [sa M il El and Ref. [H for a re- 
cent review), while Uij are the strengths of the two- 
body s-wave interactions describing atom-atom, atom- 
molecule, and molecule-molecule scattering. The diag- 
onal terms are given by Ua = Aizhaa/mi, where an 
and floo are the atom-atom and molecule-molecule s-wave 
scattering lengths. The off-diagonal terms J7oi ~ Uio 
are given by Uqi = 27rSaoi//ioi = iirhaiQ/mi, where 
aoi is the atom-molecule s-wave scattering length and 
/iQi — momi/(mo -|- mi) = 2mi/3 is the reduced mass. 

The Hamiltonian with delta-function interactions 
implicitly assumes a momentum cutoff /cmax ^ |64| . 
and we associate the coupling terms with the observed 
scattering lengths. In the numerical simulations on a fi- 
nite computational lattice, where there is always a maxi- 
mum cutoff, this description is valid as long as the lattice 
spacing in coordinate space is larger than ^ . Otherwise, 
a more careful renormalization procedure j63| is needed. 

Starting from a stable {Eq < 2Ei) molecular BEC, 
dissociation into the constituent atoms may be achieved 
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by a rapid Feshbach sweep into the atomic side of the 
resonance or by a coherent Raman transition |^ |^ 0| . 
If the transition has a negative detuning 2hA < after 
the sweep (£^0 > 2i?i), the molecules become unstable 
against (spontaneous) dissociation into free atom pairs. 
The energy level configuration after the sweep is the ini- 
tial condition for our simulations. We assume that the 
molecular BEC is initially in a coherent state, whereas 
the atoms are in the vacuum state. In addition, we as- 
sume that once the dissociation coupling x is switched 
on, the trapping potentials are simultaneously switched 
off, so that the evolution of the atomic and molecular 
fields takes place in free space. 

For molecules at rest, the excess of potential en- 
ergy 2?i|A| is converted into kinetic energy, 2?i|A| 
2?i^fc^/(2mi), of atom pairs with equal but opposite mo- 
menta around ikg, where fcp = |ko| = y^2mi\A\/h. This 
is the physical origin of the expected correlations between 
the atoms, which we will study below in greater detail 
and in the context of many-body field theory. 



To investigate the quantum dynamics we use stochas- 
tic differential equations corresponding to the positive- P 
representation of the density matrix |^ EE 1^ ■ 
The essence of this method is the mapping of the oper- 
ator equations of motion into c number stochastic differ- 
ential equations that can be solved numerically. This re- 
quires two complex stochastic fields 5'i(x, t) and $i(x, t) 
[i = 0, 1, $*(x, t) 7^ '^i{x,t)], corresponding to the oper- 
ators (x, t) and ^1 (x, t) , respectively. Ensemble aver- 
ages of the stochastic fields over a large number of trajec- 
tories (. . .)st correspond to quantum mechanical ensem- 
ble averages of normally ordered operator moments. For 
example 

([*l(x,t)]^-[*,(x',t)]") = {[M^,tt[^,i^',t)ru. (2) 

The stochastic differential equations governing the 
quantum dynamical evolution under the Hamiltonian Q 
are, in a rotating frame: 



ih 
2m 



dt 2mo 
'd7 ^ 



-V^*! - « (a + C/h$,*,) *i + X*o$i + Vx^ Ci + V-*C^oi*i*o/2 (C2 + <3) + ^J-iUll'^l d 

ih 
2mi 

VHo - i C/o.$»*.*o - 1^-? + V-*C^oi*i*o/2 (C2 - 1(3) + \/-«C/oo*o Cq, 



2mi 



■V2$o + t y C/o.$,:*.*o - + V*C^oi*i$o/2 (Ce - «C7) 

I 



iUoo^l Cio, 



(3) 



Here, the noise terms Cj iJ — 1,---,10) are real, in- 
dependent Gaussian noises with (Cj(x, t))st = and 
the following nonzero correlations (Cj(x, i)Cfc(x', <'))st = 
djkSix-x')6{t-t'). 

We note that the Hamiltonian ^ conserves the total 
number of atomic particles 



N ^ Ni{t) + 2No{t), 



(4) 



whether they are counted as bound pairs in the molec- 
ular state or as free atoms, Ni{t) = J dx>F|(x, i)>Fj(x, t) 
{i = 0,1). Since we are interested in dissociation of a pure 
molecular BEC, with no atoms present initially, then N is 
a known constant operator given by iV = 2A'o(0), where 
iVo(O) is the total initial number of molecules. Expecta- 
tion values of these operators will be denoted by removing 
the hats, e.g. N, = (N,). 



III. ANALYTIC SOLUTIONS IN THE 
UNIFORM SYSTEM 



The simplest treatment of dissociation resulting in an- 
alytic solutions is achieved under an undepleted molecu- 
lar field approximation in a uniform system T^, '13] . The 
approximation is valid for short dissociation times during 
which the total number of atoms produced is only a small 
fraction (< 10%) of the initial number of molecules. 

In this treatment we consider a uniform molecular BEC 
in a coherent state with a constant amplitude ^'o — y^po 
(where po is the uniform density), which we assume is real 
without loss of generality. The condensate is contained 
within a cubic box of side L (volume = i"^) with peri- 
odic boundary conditions, extending from —L/2 to L/2 
in each dimension. The dissociation coupling x is turned 
on suddenly, and subsequently assumed to be constant. 



The Heisenberg equations for the atomic field opera- 
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tors are then 



dt 



2mi 



9^!(x, i) . ?i „2^t A ^t ^ 

1^ ^ --i.^ — V^^'T+zAeff^'l+g*!, (5) 



2mi 



where g = X\/Po is the couphng and Aos = A + ([/oi~ 
Uoo/2)po is the effective detuning that takes into account 
the initial mean-field energy shifts due to atom-molecule 
and molecule-molecule interactions. We note that in our 
simulations we will be interested in cases with large bare 
detuning |A| 3> \Uqi — Uoo/2\po so that the mean-field 
energy shifts are negligible. 

Introducing creation and annihilation operators au_{t) 
and al^{t) for the plane- wave atomic modes with momen- 
tum k and commutation relations [ak(i)i OiL' (^)] ~ '^k,k' 
(see Appendix A), the Heisenberg equations © trans- 
form into 



dt 
dt 



2mi 



au(t) + ga] 



k' 



2mi 



Acff a 



it 



gak{t). 



(6) 



These have the following solution 



akW = Ak{t)a^{0) + Sfc(t)a^_k(0). 



where the coefficients are 



Ak{t) = cosh{gkt) - iAk smh{gkt)/gk, 
Bk{t) = g smh{gkt)/gk, 



(7) 



(8) 



and we have introduced A^ = fi,k^/(2mi) -I- Aoff and 
gk = [g^ - A^]i/2, where fc = |k|. The coefficients Akit) 
and Bk{t) satisfy \Ak{t)\'^ - Bl{t) = 1. 

Given the initial vacuum state of the atomic field we 
can calculate any expectation values of the field opera- 
tors. In particular, we find that the only nonzero second- 
order moments are the particle number per mode and the 
pairing field per mode: 

n^{t) = (nk(t)) = Bl{t) = [g/gkf sinh'igkt), (9) 



mk(t) = (mkW) =Ak{t)Bk{t), 



(10) 



where h\^{t) — a\j^t)a-\s_{t) and TOk(i) — ak(i)a-k(0 ^''^ 
the respective operators. All other second-order mo- 
ments are equal to zero. Higher-order moments factor- 
ize according to Wick's theorem and can be expressed 
in terms of the above second-order moments. From 
|^fe(OP ~ B\{t) = 1, we find in addition that n\^{t) and 
mk(i) are related by 



|mk(t)P =nk(t)[l+nk(t)]. 



(11) 



Using these solutions, the total number of atoms pro- 
duced is 



N,{t)=Y^ n^{t)=Y^ Bl{t). 



(12) 



The treatment of an infinite system in free space 
is achieved as usual by taking the limit of L — > oo 
(Afc 0) and transforming the mode creation and an- 
nihilation operators into continuous Fourier transforms 
a(k) and a^(k) that satisfy the commutation relation 
[a(k), a^(k')] = 6{k — k') (see appendix A). In this case, 
the solutions for the normal and anomalous moments are 



(a^(k,t)a(k',i)) 
(a(k,t)a(k',i)) 



Bl{t)5{k-k'), 
Ak{t)Bk{t)5{k- 



k!) 



(13) 
(14) 



Atom-atom correlations 



We are interested in the atomic correlations result- 
ing from momentum conservation in the molecular dis- 
sociation. These may be quantified using a number of 
different, but related, density-density correlation func- 
tions. Furthermore, these measures may correlate fields 
at points in position space or momentum space. Since the 
fundamental correlation is between momenta, the mo- 
mentum correlation functions are simpler to interpret. 
However, experimentally it is easier to measure spatial 
correlations after time-of-flight expansion of the cloud. 
In the present paper we will concentrate on analyzing 
atomic correlations in momentum space. Spatial corre- 
lations after expansion will be studied in a future work 

Starting from the atomic vacuum, and assuming no ini- 
tial momentum spread in the molecular BEG and no scat- 
tering, molecular dissociation produces pairs of atoms 
with equal and opposite momenta. In the short time limit 
this is well approximated by the analytic solutions of the 
previous subsection. The strength of the atomic correla- 
tions may be quantified via Glauber's second-order cor- 
relation function [491 



ff(2)(k,k',t) = 



(4W4'W«k'(^)ak(^)) 



(15) 



It is defined in terms of normally ordered operator prod- 
ucts and is normalized so that it is dimensionless. The 
pair correlation describes the ratio of the probability of 
joint detection of pairs of atoms with k and k' to the 
product of probabilities of independent atom detection 
events at k and k'. For example, i?'-^-' = 1 for uncorre- 
lated atoms, g^^-* = 2 for thermally bunched atoms, and 
g(2) > 2 indicates super-thermal bunching. 

Note that the conditions for an expansion using Wick's 
theorem are met by the solutions for the uniform sys- 
tem under the undepleted molecular field approxima- 
tion, Eq. ©. In this case, the pair correlation function 
(7^^^ (k, k', f) can be factorized and be expressed in terms 
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of products of the second-order normal and anomalous 
moments (ajj.(i)ak' (t)) and (ak(i)ak' (i))- For example, 
in the case of k' = — k, the pair correlation function is 



7(2)(k,-k,i) = 1 



1 



= 2 



(16) 



where we have taken into account that (nk(i)) = 
(n_k(t)), {al^{t)a-k{t)) = and used Eq. For other 

pairs of momenta the correlation function is obtained in 
a similar way, with the overall result given by 



1, k'^k,k'^ 

2, k' = k 7^ 0, 

2 + l/nk(t), k' = -k, 

3 + l/nk(i), k' = k = 



(17) 



As we see, the correlation function for atom pairs with 
equal but opposite momenta is the strongest [except for 
the special case of k = k' = 0]. Quantitatively, strong 
superbunching in (^^^^(k, — k, t) is achieved for low mode 
populations, nk(t) <C 1. In this regime, the numerator 
in Eq. (|15|l is approximately proportional to the mode 
population n\^{t), while the denominator is the product 
of populations, and therefore g'^^^(k, — k, t) « l/nk(t) '3> 
1. As the mode occupancies grow with time, the pair 
correlation approaches the level of thermal bunching 

Thus, at high densities Glauber pair correlation be- 
comes a less sensitive measure of the correlation strength. 
This is despite the fact that the correlation between the 
atoms with equal and opposite momenta is still maximal 
at any given density, in this analytically soluble model. 
There is perfect squeezing of particle number-difference 
fluctuations below the shot-noise level. This is quantified 
via the normalized variance 



F(k,k',i) = 



(nk(i)> + (nk'(t)) 



[A (nk(^)-nk'(^))] 
(nk(<)) + (nk'(i)) 



(18) 



where AA ^ A — (A) is the fluctuation and the colons 
:: indicate normal ordering of the creation operators be- 
fore the annihilation operators. This definition uses the 
conventional normalization with respect to the shot-noise 
level resulting from a Poissonian number probability dis- 
tribution, such as for a coherent state. This is given 
by the sum of the mean occupation numbers, {nu{t)) + 
(nk'(i)). Variance smaller than one, F(k, k',t) < 1, im- 
plies reduction of fluctuations below the shot-noise level 
and is due to correlation between particle number fluc- 
tuations in the k and k'-modes, while F(k, k',i) = 1 for 
uncorrelated modes. 

For the solutions of Eq. ® the number-difference 
variance for equal but opposite momenta is given by 



V{\s., — k, <) = at all times, i.e., for all occupation num- 
bers nk(i) — ri_k(i). This is the ideal case, implying per- 
fect correlation between fluctuations in rik(t) and n_k(t) 
and corresponding to 100% squeezing below the shot- 
noise level. 

We note that for two equally occupied modes, the vari- 
ance V^(k, — k, t) and the second-order correlation fimc- 
tion g^^^ (k, — k, t) are related by 

V{k,-k,t) = l~nk(i)[.9^'Hk,-k,t)-g(2)(k,k,t)], (19) 

where we have used g*^^^(k, k, t) — (7^^'(— k, — k, i) due 
to the spherical symmetry of the problem. Therefore, 
perfect noise reduction, V{k,—k,t) = 0, implies that 
the maximum degree of correlation is gn^ax{k,—k,t) = 
g^'^\k,k,t) + l/nk(i). Using the thermal level of auto- 
correlation that occurs in the analytic solutions, Ea. (|17|l . 
g'-^\k,k,t) = 2, this gives 



ff^?L(k,-k,0 = 2 + l/nk(t), 



(20) 



which coincides exactly with the analytic correlation, Eq. 

m- 

In a similar way, one can show that the maximum 
(ideal) degree of correlation Eq. (|2U|I corresponds to the 
maximum anomalous moment mk(i) 



max{|mk(t)P} = nu{t)[l + nu{t)] 



(21) 



which again coincides with the solution, Eq. Hll|l . of 
the present analytically soluble model. This result fol- 
lows immediately from the Wick-factorized expression for 
5(2)(k,-k,t), Eq. 



j'^^Hk,~k,t) ^ 1 



(22) 



assuming (dj^(i)a_k(t)) = 0. 

More generally, i.e., for nonideal cases, the number- 
difference variance V{k, — k, t) would be larger than zero 
and can serve as a sensitive measure of the correlation 
strength between particle number fluctuations in differ- 
ent modes. In the nonideal cases, V{k,—k,t) > also 
implies that the following inequalities hold: 



ff(2)(k,-k,i) < 2 + l/nk(i), 

|mk(i)|' < nkit)[l + nk{t)l 



(23) 
(24) 



provided that one still has ^'^^^(k, k, t) = 2 and 

(4(t)a_k(0) = 0. 

Quantifying the strength of correlations via the 
number-difference variance V(k,—k,t) is especially use- 
ful at large mode occupancies nk(i) 3> 1, when the 
pair correlation approaches (;^^^(k, —k,t) — > 2 and there- 
fore becomes a less sensitive measure of the correlation 
strength. 
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IV. NUMERICAL RESULTS AND DISCUSSION (a(k, t)a(-k, t)) [having units of m^] via 



Our goal in the analysis of atom correlations in molecu- 
lar dissociation is to go beyond the undepleted molecular 
field approximation and to treat the quantum dynamics 
for realistic nonuniform systems. Analytic solutions for 
these cases are no longer available, and we will resort to 
numerical solution of the positive-P equations (j2Jl, which 
describe all quantum effects from first principles. 

Using the positive-P method, we are able to extend the 
analytic results of the previous section in three ways, to 
include: (a) the depletion of the molecular condensate; 
(b) s-wave scattering interactions; and (c) nonuniform 
condensates. 

In the first two cases, we will consider uniform systems, 
while the nonuniform case will be simulated within the 
undepleted molecular field approximation. The reason 
for doing this is to give a precise quantitative descrip- 
tion of each of these effects and to assess their relative 
importance in altering the atom-atom correlations. In 
addition, the results obtained in each case can serve as 
a benchmark for other approximate numerical methods. 
An example is the pairing mean-field method of Ref. 
which treats the molecular field depletion at the level of 
a coherent state, while the dynamics of the atomic field 
is treated via the normal and anomalous densities. The 
range of validity of this method will be examined below 
using the exact results of the present treatment. 

The positive-P numerical method is implemented on 
a uniform spatial lattice with periodic boundary condi- 
tions, and the continuous fields are therefore represented 
by discrete mode amplitudes. The momentum space lat- 
tice is reciprocal to the spatial lattice, with Ak = 2'k/L 
being the lattice spacing and L the length of the spatial 
domain in each dimension. 

Since the numerical method is able to treat both uni- 
form systems in a finite box as well as nonuniform sys- 
tems in free space, the momentum-space field amplitudes 
are commonly treated via the lattice-discretized momen- 
tum components a(k, t) and a^(k, i), which correspond 
to the continuous Fourier transforms (see Appendix A) 
in the limit Afc ^ (P ^ oo). (The correspondences 
between the operator Fourier components and those of 
the stochastic fields are outlined in Appendix B). Ac- 
cordingly, for any finite computational lattice one can 
formally identify a set of momentum modes described 
by creation and annihilation operators a\^{t) and a\j^t) 

(with commutators [a]^{t),a\.,{t)\ = (Sk,k')' which are 
related to their lattice-discretized continuous counter- 
parts a(k, i) and a^(k, i) via d\^{t) = a(k, i) (Afc)"^^^ and 
a],(t) = at(k,t) {Akf\ 

Similarly, the (dimensionless) mode populations nk(t) 
and the pairing fields per mode TOk(t) are related 
to the normal and anomalous densities n(k, = 
(7T,(k,i)) = (at(k,t)a(k,t)) and ■m{^^,t) = {■m{^^,t)) = 



nk(t) = n(k,<) {Ak) 
mi^(t) = m(k,i) (Afc)^ 



(25) 
(26) 



In the continuous limit, these should be understood as 
corresponding to the expectation values of the respective 
operators defined as integrals over the momentum-space 
lattice volume element w(k) — (Afc)"^ around k, 

hi^{t)^[ h{k',t)dk' ~h{k,t){Ak f , (27) 
rhkit) = / rh{k', t)dk' ~ m(k, t) {Ak f , (28) 

JD(k) 

where Afc is small enough so that n{k',t) and 7fik(k',t) 
under the integrals do not vary much within the integra- 
tion volume and can be replaced by ?T.(k, t) and m{k, t). 
The total number of atoms is 



Ni{r) 



J n(k,t)dk ~^^n(k,t) (Afc)' 



(29) 



The normalized second-order correlation function is 
now defined via 



g^'\ky,t) = 



{k,t)aHk' ,t)a{k' ,t)a{k,t)) 
(n(k,<))(n(k',t)) 
(:n(k,<)n(k',t) :) 



(n(k,<))(n(k',t)) ' 



(30) 



which is equivalent to Eq. H15|l . provided that a(k, t) and 
ak(i) are defined on the same momentum-space lattice. 

Finally, we define the normalized variance of the par- 
ticle number difference rik(t) — ?^k'(^)J 



K(k,k',i) = 1 



[A (nk(0-nk'(^))]' 
(nk(t)> + (nk'(t)) 



(31) 



where the particle number operators are to be understood 
according to Eq. (|27|l . The subscript v in Eq. H31|l signi- 
fies the fact that in the continuous model describing an 
infinite system, the variance is defined for the counting 
volume V — (Afc)'^, and therefore Eq. is equiva- 

lent to Eq. ((TH)l . Accordingly, the results for Vv{k,k\t) 
depend on this volume; unlike the pair correlation func- 
tion (7^^^ (k, k', f), the scaling of the numerator and of 
the denominator with respect to {Ak)^ does not lead to 
cancellation. The explicit dependence on (Afc)"^ becomes 
evident if we rewrite Eq. I|31(l via the densities n(k, t), 
which are independent of (Afc) : 



K(k,k',i) = l + (Afc)' 



[A{n{k,t) -n{k\t))]' 
(n(k,t))-|-(n(k',t)) 



(32) 

As we see, for a given negative value of the normally 
ordered variance (: [A(n(k, i) — n(k',t))]^ :), the degree 
of squeezing below the shot noise level degrades {Vy be- 
comes closer to one) as the counting volume is decreased. 
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Specific numerical results obtained for 14(k, k',t) 
should not be confused with the requirement that phys- 
ical observables should not depend on the choice of the 
lattice spacing Afc. Rather, the calculated variance for a 
given computational grid corresponds physically to fluc- 
tuations in the difference of particle numbers in the lat- 
tice volume element (Afc)'^. For a given density, the num- 
ber of particles in a smaller volume is lower and the corre- 
lations between number-difference fluctuations becomes 
weaker. Stronger correlation signals can be obtained via 
the procedure of binning (see below), in which case the 
counting volume is chosen to be larger than (Afc)^. 

To further clarify this point we note that the unnormal- 
ized variance of the continuous iie?7.sit2/-difference fluctu- 
ations, for k 7^ k', 

([A (n(k, t) ~ n(k', t))f) = [(n(k, t)) + (n(k', i))],^^^) (Q) 
+ (:[A(n(k,<)-n(k',t))]':), (33) 

includes (after rewriting it in terms of the normally or- 
dered operator products) a term proportional to a delta- 
function ^'^^^(0), which is the shot-noise level. While 
strong quantum correlation and the reduction of density- 
difference fluctuations below the shot-noise level is gener- 
ically due to negative values of the normally ordered vari- 
ance, the degree of squeezing below the delta-function 
shot-noise is not well defined via Eq. 133|) . On the other 
hand, operating with the variance of fluctuations between 
the particle numbers, as in Eq. H31I) , makes the definition 
of the shot-noise unambiguous and therefore the degree 
of squeezing is well defined. Albeit, the squeezing now 
depends on the counting volume (Afc)'^. 

A. Parameter values 

Taking a ^''Rb (mi = 1.43 x 10"^^ kg, toq = 2mi) ex- 
periment as our specific example, we have performed sim- 
ulations for the following set of parameters. The initial 
molecular BEG density was chosen to be po = 5 x 10^^ 
m~^ for the uniform system. For the nonuniform sys- 
tems, the same value was chosen to be the molecular BEG 
peak density. The size of the uniform system L was the 
same in all three spatial dimensions, with L = 1.38 x 10~^ 
m, and the lattice grid contained 64'^ points. Thus, the 
lattice spacing in momentum space was Afc = 271 /L = 
4.56 X 10^ m^^, while the maximum cutoff momentum 
was fcmax = 1-46 X 10^ m^^, which is about twice the res- 
onant momentum fco = 1/2777,1 |Acff|/?i and exceeds the 
physically relevant range of momenta. The total initial 
number of molecules was A^o(O) — poL^ — 1.3 x 10^ in 
all uniform cases. (The values for the nonuniform system 
are given m Sec. CSdJ) 

The atom-molecule coupling strength was x = 7 x 10^^ 
m^/^/s. The bare detuning A in different cases was 
chosen to result in the same initial effective detuning 
Aoff = A-|-([/oi-C^oo/2)po = -1.96x10'^ s-^ irrespective 
of the mean-field interaction contributions. The values 




FIG. 1: Fractional population in the molecular and atomic 
fields, No{t)/No{0) and Ni{t)/2No{0), as a function of time 
t. The time is in units of to = I/5 = l/(X\/Po)' '^hich in this 
example is to = 0.2 ms. The dash-dotted line is the fraction 
of molecules in the noncondensate modes, k 7^ 0. The dashed 
line is the analytic result for the total atom number in the 
undepleted molecular field approximation. 

of s-wave scattering lengths for the molecule-molecule, 
atom-molecule, and atom-atom interactions are given in 
the respective subsection below. In all simulations, the 
number of stochastic trajectories for calculating the ex- 
pectation values was 7500 (unless stated otherwise) and 
the time step was At = 1.26 x 10~^ s. Most of the sim- 
ulations are performed for durations t/to = 2.5, where 
to = {x^/po)~^ = 0.2 ms is the time scale. 

With a spatially uniform molecular field and no de- 
pletion or two-body interactions, we validated our code 
by comparing the numerical results with those of the an- 
alytically soluble model, Eq. jS)). The numerical code 
for solving the positive-P stochastic differential equations 
was implemented using the XMDS software package f65l| . 
Within the sampling errors, the simulated momentum 
space mode populations 7ik(t) — r7(k, t) (Afc)'^ and the 
correlation functions were in excellent agreement with 
the analytic results given by Eqs. lO, H10|) and H17|l . 
The subsections below extend these results to include the 
molecular depletion, s-wave scattering interactions, and 
the treatment of nonuniform systems. 

B. Role of the molecular field depletion 

1. Total particle numbers and atom-atom recombination 

Here, we simulate Eqs. Q for a uniform system and 
all s-wave scattering interactions set to zero (ay = 0, 
and hence A — Acff). Including the molecular field in 
the dynamics leads to its depletion and a corresponding 
reduction in the number of atoms produced compared to 
the analytically soluble case. 

In Fig. n we plot the total fractional number of 
molecules and atoms, No{t)/No{0) and 7Vi(t)/27Vo(0), as 
a function of time. The dashed line shows the total atom 
number in the undepleted molecular field approximation. 



Eq. H12() . for comparison. We see that the undepleted 
molecular-field result agrees well with the exact quan- 
tum result for time durations resulting in the conver- 
sion of less than 10% molecules. In the present example, 
10% conversion corresponds to 2.6 x 10^ atoms produced 
[Ni{t)/2No{0) ~ 0.1] and occurs at t/to ~ 1. At this time 
the discrepancy between the exact numerical and the ap- 
proximate analytic results is ~ 9%. For longer times 
the discrepancy increases and the approximate analytic 
result eventually produces an exponentially growing out- 
put. This unphysical behavior is an obvious consequence 
of the fact that the undepleted molecular field approxi- 
mation is no longer valid. 

Physically, the total number of atoms produced satu- 
rates due to the depletion of the molecular condensate 
containing a finite number of molecules to start with. At 
the same time, the dynamics of dissociation - after the 
initial spontaneous regime - is affected by the reverse 
process of recombination of atom pairs into molecules. 
The dashed-dotted line in Fig. Q shows the total num- 
ber of molecules outside the condensate mode k = 0. 
Since the dissociation in this uniform system starts with 
all molecules being initially in the condensate mode, the 
population of the k 7^ modes can only occur (in the 
absence of elastic collisions, Uij = 0) due to atom-atom 
recombination. Note that the recombination can in prin- 
ciple involve atom pairs with any momenta, and not nec- 
essarily those with equal and opposite momenta. This 
has the effect of reducing atom-atom correlations in the 
long time limit (see below), hence the term "rogue asso- 
ciation" which we will use. 

The population of the molecular noncondensate modes 
implies effective heating of the molecular gas and can be 
regarded as the reverse of the rogue dissociation known to 
occur in the opposite process of association of an atomic 
BEC into a molecular BEC 62, 66] . In the example of 
Fig. n 18% of the remaining molecules at t/to ~ 2.5 
are outside the condensate mode. At the same time, the 
total number of molecules remaining is about 45% of the 
initial number. 



2. Atomic momentum distribution 

In Fig. |2Ja) we plot a slice through the origin of the 
atomic momentum density distribution ?i(k, t) at t/to = 
1. We see a clear ring structure corresponding to the 
distribution of atoms around the surface of a sphere of 
radius feg. Similar ring structures have been observed 
in the time-of-flight spatial column densities of atoms 
produced in dissociation experiments using ^^Rb2 and 
molecules 013. 

To simplify the presentation and comparison of the 
results we make use of the spherical symmetry of 
the problem and also plot angle-averaged distributions 
(n(k, t))^^g and {n\^{t))^^g [see Fig. EJb)] as a function of 
the absolute momentum ]k]. Here, the brackets (. . 
refer to the procedure of averaging of the quantum me- 
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(a) n(k,y [^im^] 




FIG. 2: (a) (color online) Slice through the origin of the 3D 
atomic density in momentum space 7i(k, t) [in units of /im^] 
a.t t = to, where to — l/(X\/Po) ~ ^■■^ momentum 
components kx and ky are in units of ko = y'2mi|Acfr|/ft, 
which in this example is ko = 7.3 x 10® m~ . (b) Angle- 
averaged density distribution {n(k, t))^,e and mode occupa- 
tions {nk(t))ip,g = {n(k, i))^,e(Afc)^ as a function of the abso- 
lute momentum |k| at t = to (solid line). The lattice spacing 
is Ak — 4.55 X 10^ m~^. The dashed line is the result for the 
analytically soluble case with the undepleted molecular field, 
Eq. @. 



chanical expectation values over the angles ip and in a 
spherical coordinate system. 

Comparison of the exact numerical result for the angle- 
averaged population distribution (nk(t))i^.e and the re- 
spective analytic result, Eq. 0, shows that the discrep- 
ancy at t/to = 1 is less than 8%, with the peak value of 
{n\i{t))^fi at ]k]/fco = 1 being reduced from the unde- 
pleted result of sinh^(l) — 1.38 to 1.27. For longer time 
durations the discrepancy increases. The small oscilla- 
tory structure seen in the wings of the distribution func- 
tions (for both the exact numerical and analytic curves) 
is additional evidence that the sampling error due to 
stochastic averaging is small. This oscillatory behavior 
is characteristic of the highly detuned modes outside the 
spherical shell as these modes do not experience an ex- 
ponential growth but rather oscillate at the spontaneous 
level. In terms of the approximate analytic solutions, the 
oscillations take place for the modes with k — ]k] for 
which the gain coefficient gu = [g^ — A^]^/^ is pure imag- 
inary. This occurs when = [h]s? / (2mi) -I- Acfi]^ > 5^, 
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FIG. 3: (a) (color online) Slice through the origin of the 3D 
correlation function g^-^\k., —k.,t) ai t — to. The correspond- 
ing density profile is shown in Fig. HJa) . The sampling er- 
rors due to stochastic averaging are small only within the 
region of 0.65 < |k|/fco < 1.25, where the correlation func- 
tion ranges between ~ 2.8 and 6. (b) Angle-averaged cor- 
relation functions {(j'^' (k, k', (for k' — — k, k, kgo) at 
t = to as a, function of the absolute momentum |k| within 
the spherical shell 0.65 < |k|/fco < 1.25. The solid line with 
the dots (right scale) is the respective angle-averaged variance 
{K(k,-k,t))^,9. 



in which case the sinh^(gfci) = sinh^(z|(7fc|t) term in Eq. 
turns into sin^(|gfe|i) thus producing the oscillations. 



3. Atom-atom correlations 

The pair correlation function for the atoms with op- 
posite momenta, g*^^^ (k, — k, t), is shown in Fig. Ola). 
Here, we plot a slice through the origin of the 3D cor- 
relation function at t/to — 1. Due to the finite num- 
ber of stochastic trajectories and the normalization of 
the correlation function to the product of mode occupan- 
cies, the stochastic ensemble average gives highly noisy 
results in the regions of kspace where the mode popula- 
tions are small (nt ^ 0.3). This is seen within the central 
(|k|/A;o < 0.65) and outer (|k|/fco > 1-25) region of the 
slice plot where the sampling error is large |67l |. Within 
the spherical shell of 0.65 < |k|/fco < 1.25, on the other 
hand, the mode populations are higher and the sampling 
errors are negligible. Within this region the results for 



g^^^ (k, — k, t) scale according to the colormap shown on 
the right of Fig. E^a). 

Figure|2Ib) shows the angle-averaged correlation func- 
tions (g'^^^ (k, k', t));p.e within the spherical shell 0.65 < 
|k|/A;o < 1.25, at t/to — 1- As expected, for k' = — k the 
correlations are super bunched, {g^-^\k,—'k,t))^^0 > 2, 
while for k' = k we see the thermal level of bunch- 



2 characteristic of Gaussian 



ing (5(2)(k,k,t))^,e 
statistics. The line showing an uncorrelated level of 
(g^^-* (k, kgo, i)),^,e = 1 is for the pair of momenta k and 
k' = kgo , where kgo corresponds to a 90° degree rotation 
of the density distribution in the {k^ky) plane around 
the kz axis. 

Apart for the fact that the pump depletion results in 
slightly decreased mode populations relative to the ana- 
lytic result, the correlations follow closely Eqs. H17|l . ex- 
cept that nk(i) is now the actual, decreased, mode popu- 
lation. This observation is valid for at least the simulated 
time window of t/to — 2.5, corresponding to about 55% 
conversion. This is expected since the correlations are 
due to momentum conservation, which is unaffected by 
the decreasing rate of conversion, which is the major ef- 
fect of the molecular depletion. Thus, we find that the 
pair correlation for equal and opposite momenta is still 
well approximated by 



(g(2)(k,-k,t))^,(,^2 + l/(nkW)^,e. 



(34) 



Since {n]s_(t))^_g is slightly less than in the unde- 
pleted molecular field approximation, the pair correlation 
Eq. (|34|l is accordingly higher. In Fig. Otb) the result of 
Eq. (|34|l is plotted by the dotted line and is almost in- 
visible as it follows closely the actual numerical value 
calculated using Eq. (|Sn)l . 

The fact that the correlations can be approximated 
by Eq. (|34|l shows that Wick factorization continues to 
approximately hold for time durations corresponding to 
~ 55% conversion. We recall that the ingredients of this 
result are the vanishing off-diagonal normal moments, 
{a\^(t)(X-]s_{t)) = 0, and the relationship of Eq. lfTT|l be- 
tween the normal and anomalous moments. These ingre- 
dients continue to hold in the present example. 

Figure ^ shows the temporal behavior of the angle- 
averaged pair correlation functions (t;'-^-' (k, k', i))^^e at 
|k|/fco = 1 corresponding to the peak momentum in 
the density distribution. The approximate result for 
(.g^^^(k, — k, i))y^6(, given by Eq. is plotted with a 

dotted line, but is coincident with the actual numerical 
value shown by the solid line. 



4- Number- 1 



Terence variance 



As we mentioned earlier, the pair correlation func- 
tion becomes a less sensitive measure of the correlation 
strength when the mode occupation numbers increase. 
In particular, as nk(t) increases with time the pair cor- 
relation approaches (7*^^-'(k, — k, t) 2 which coincides 
with the level of thermal bunching for the autocorrelation 
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FIG. 4: Dynamics of the angle- averaged pair correlation func- 
tions ((;'^'(k, k', at |k|/fco — 1, for k' = — k, k, and 
kgo. The time is in the same units as in Figs. Q and |21 
The solid line with the dots (right scale) is the variance of 
the number-difference fluctuations for equal but opposite mo- 
menta, (K(k, — k, at |k|/fco = 1. 



function. In this regime, the variation in ^(^^(k, — k, t) is 
no longer significant. As a result, distinguishing between 
strong and weak correlations becomes problematic. 

A more sensitive measure of the correlation strength 
is provided by the variance of the particle number- 
difference fluctuations between the correlated modes, 
Eq. H18|l . The angle- averaged number-difference vari- 
ance for the modes with equal but opposite momenta, 
(K(k, — k, is shown in Figs. Elb) and ^ by 
solid lines with dots. As we see from Fig. ^h), the 
variance is suppressed well below the shot-noise level 
(K(k, —k,t))tp^9 < 1 for the entire spectral range. 

Similarly, the variance {Vv{k, —'k,t))^^g shows sub- 
shot-noise fluctuations in the entire simulated time win- 
dow (see Fig. It increases above the ideal result of 
(14(k, — k, — towards the end of the simulated 

time window. In particular, at t/to = 2.5 it is given 
by {Vy{'k,—'k,t))^^g ~ 0.22 which corresponds to 78% 
squeezing below the shot-noise level. The reduction of 
squeezing is attributed to atom-atom recombination into 
molecules with nonzero momenta, k 7^ (see Fig. 
Since this process involves pairs of atoms with nonoppo- 
site momenta, a certain fraction of the atoms is left with- 
out their correlated partners of opposite momentum. As 
a result, the overall strength of the correlation signal for 
the atoms with equal and opposite momenta is reduced. 

The reduction of the correlation strength due to atom- 
atom recombination becomes a sizeable effect after dis- 
sociation times {t/to > 2) corresponding to more than 
^ 35% conversion. At very short times (less than 10% 
conversion) our exact numerical results confirm that the 
predictions of the analytically soluble model with the un- 
depleted molecular field are a good approximation. At 
longer times the correct treatment of the system requires 
the inclusion of the molecular field depletion. 

As shown in Ref. the molecular depletion can be 
taken into account using a pairing mean-field method 
which is simpler to implement numerically than the 
present (exact) positive-P method. In the mean-field 



aij (nm) 


Uij (xlO"^'^ m^/s) 


tsim/to {to = 0.2 ms) 


aoo = 2 


Uoo ~ 0.921 


2.125 


aoi = 2 


Uoi = 1.38 


1.75 


Oil = 2 


Uu = 1.84 


1.25 


ci^j — 2 


Uij - as above 


1.25 


Oil = 5.3 


Uii = 4.88 


0.5 



TABLE I: Scattering lengths aij, coupling constants Uij, 
and the respective simulated time duration isim (in units of 
to — l/(X\/Po) ~ 0.2 ms) for simulations with the s-wave 
scattering interaction terms. For each case, the values of the 
remaining terms were set to zero, i.e., aij = except for the 
quoted term. All other parameters are as in Sec. IV A. 



method, the molecules are described as a coherent field 
in the condensate mode (&k=o(0) while the atoms are 
treated at the level of diagonal normal and anomalous 
populations, (aj^(t)ak(t)) and (ak(i)a-k(i))- While im- 
posing Wick's factorization scheme and assuming perfect 
correlation between the atoms with opposite momenta at 
all times, this treatment does not allow for atom-atom re- 
combination into molecules outside the condensate mode. 
Accordingly, the predictions of the theory are a good ap- 
proximation only for dissociation durations resulting in 
no more than ^ 35% conversion, as can be seen from the 
comparison with the present exact results. 

The pairing mean-field method performs better in the 
long time limit if the only observables of interest are the 
total particle numbers. Its performance can be expected 
to improve for calculating the correlation functions if 
it is expanded to include the molecular noncondensate 
modes and the off-diagonal terms in the atomic normal 
and anomalous densities. 

Before turning to the analysis of atom correlations in 
the presence of s-wave scattering interactions, we note 
that the positive-P simulations without the two-body 
collisions fail for durations longer than t/ta ~ 2.5. At this 
time about 55% of molecules have been converted into 
atoms. The failure of the simulations for longer durations 
is due to the generic boun dary term problem that occurs 
in the positive-P method [53,|53- Despite this, our re- 
sults are a significant improvement and extend beyond 
the regime of validity of the simple undepleted molecular 
field approximation or the pairing mean-field theory. A 
possible route to increase the simulation time is to use 
stochastic gauges jl^ls^, which is, however, beyond the 
scope of the present paper. 



C. Role of the s-wave scattering interactions 

Here, in addition to simulating the molecular field 
depletion we include the s-wave scattering interaction 
terms. Table shows the different cases considered, to- 
gether with the simulated time window tsim- The sim- 
ulations are shorter than in the absence of the s-wave 
scattering since the positive-P boundary term problem 
|55l Is^ becomes more severe for quartic interactions. 
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The simulation times were chosen to be less than those 
for which the sampling error was observed to grow rapidly 

The shortest tsim is for the case of an = 5.3 nm 
which is the natural background scattering length for 
^^Rb atoms. Longer tgim are achieved with smaller values 
of aij. The scattering length can in principle be tuned 
during dissociation by applying a magnetic field in the 
vicinity of a Feshbach resonance. In this case the dis- 
sociation itself would have to be invoked independently, 
using optical Raman transitions. While this is unlikely to 
simultaneously result in the specific values of aij quoted 
in Table m our main goal is the understanding of the role 
of each physical process separately - by quantitatively 
estimating their relative importance and the potential 
disruptive effect on atom-atom correlations. 

The results of simulations of the cases of Table m show 
that the momentum distribution and the pair correla- 
tion for atoms with opposite momenta remain (within 
the simulated time windows) very close to the values ob- 
tained with no s-wave interactions. For example, in the 
case of On — 2 nm, the peak occupancy reduces from 
1.70 to 1.62 at t/to = 1.25, while the corresponding pair 
correlation still follows Eq. (|34|l . The reduction in the 
mode population is attributed to the mean-field energy 
shift due to atom-atom interactions, which dynamically 
changes the phase matching condition and hence results 
in a drift of the effective detuning Acg as the populations 
evolve. As a result, the resonant momentum also drifts 
causing less efficient conversion into the initially resonant 
mode. Nevertheless, since the bare detuning |A| in our 
examples is still much larger than all mean-field phase 
shifts Uij{^j{x,t)^j{x.,t)), the resonant momentum is 

still well approximated by fep — y^2mi\A\/h. 

The expected reduction of the actual correlation 
strength due to s-wave collisions and the resulting re- 
distribution of the atomic momenta is best revealed 
through the variance of the number-difference fluctu- 
ations, {Vu{]<., —'k,t))^g. For example, in the case of 
aij = 2 nm (which is the worst case, among the sim- 
ulations surviving up to durations of t/to = 1.25), the 
variance at the resonant momentum fcg increases from 
the colhsionless resuh of 0.024 to 0.032 at t/to = 1-25. 
This implies the reduction of number-difference squeez- 
ing from 97.6% to 96.8%. The overall effect is, however, 
smaller than decorrelation due to atom-atom recombina- 
tion discussed earlier, which in this example is responsible 
for the reduction of squeezing from 100% down to 97.6%. 

D. Nonuniform systems and mode mixing 

We now turn to the analysis of atom-atom correlations 
in nonuniform systems corresponding to dissociation of 
realistic trapped BECs. The net effect of inhomogeneity 
is mode mixing, which can dramatically affect the cor- 
relations. The quantitative details are best understood 
in the undepleted molecular field approximation. In this 



case, the atom-atom recombination and s-wave scattering 
interactions are absent, and any reduction in atom-atom 
correlation strength - compared to the uniform case - is 
due to mode mixing. 

Hence, the molecular field amplitude ^'o(x, 0) = 
y'/DQ^x) [where jOo(x) is the density] can be absorbed into 
an effective coupling .g(x) = x\/ Po(x). This leads to the 
same equations of motion for the atomic field operator 
\['i(x, i) as Eq. except that the coupling constant 
g is now a function of x. Introducing the Fourier trans- 
form ^(k) = (27r)^'^/^ / c?xg(x) exp(— ik-x), the equation 
of motion for the atomic Fourier component a(k, t) is 

daCk,t) fm^ . \.,, , 

+ (2^ / dk'~gik + k')a\k',t). (35) 

As we see, a(k, t) couples to the range of momentum 
components k' in the conjugate field a''{k',t), which 
we refer to as mode mixing. This is in contrast to 
the uniform case, where ^(k + k') is the delta function 
g{k + k') = (27r)3/2gJ(k + k'), so that a{k,t) couples 
only to the conjugate field at the opposite momentum 
a^-k,t). 

The overall effect of the mode mixing is to correlate 
the atoms with momentum k not only with the opposite 
momentum — k, but also with the atoms distributed in a 
range of momenta around — k. As a result, for a given 
strength of g{k) the pair correlation between k and — k 
is expected to be reduced compared to the uniform case. 

To characterize mode mixing quantitatively we simu- 
late the positive-P stochastic differential equations that 
are equivalent to Eqs. H35|) . The initial molecular BEC 
is assumed to be formed in a spherically symmetric 
harmonic trap with equal trap oscillation frequencies 
uj = LiJx,y,z- It is assumed to be in a coherent state, 
with the density distribution ^iven by the Thomas-Fermi 
parabola po(x) = po(0)(l— |x| /Rxf) for |x| < Rtf, and 
po(x) — for |x| >Rtf- Here, po{0) is the peak density, 
while Rtf — \/2fi,[/ooPo(0)/(TOoW'^) is the Thomas- Fermi 
(TF) radius. 

1. Case 1: weak inhomogeneity 

The choice of parameter values in this set of simula- 
tions is specifically targeted to give a situation which 
is directly comparable with the results of the previ- 
ous uniform system. The uniform results for a cubic 
box of side can be regarded as an approximation 
to a realistic nonuniform system if L„ is matched with 
the characteristic size ~ 2Rtf of the initial molecu- 
lar BEC, while the uniform density po is matched with 
the peak density po(0). More specifically, we choose 
Rtf = (877/15)^^/"^ L„ and po(0) — po, which gives 
the same initial total number of molecules in both cases, 
No{0) = (87r/15)po(0)i?3p = poLl 
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FIG. 5: Comparison between the results of simulations of weakly inhomogeneous {Rtf = 11.6 fim) and strongly inhomogeneous 
{Rtf = 4.82 fim) systems, columns 1 and 2, having the same peak density po(0) = 5 x 10^® m~"^. Column 3 refers to the 
same physical system as in column 2, except that the length of the simulated spatial domain L is halved. Note that the lines 
are provided as guides-to-the-eye between calculated data points, (a) Slice through the origin of the initial molecular BEC 
density profile, po(2;,0, 0). (b) Angle -averaged momentum density (?2(k, and distribution of mode occupations {Tiu.(t^)^p^6 

at t = to as a function of the absolute momentum |k|/fco, where ko = 7.3 x lO" m~^ as in Fig. |5| (c) Fraction of the total 
number of dissociated atoms relative to 2A'^o(0) vs time, (d) Angle- averaged mode population {n(k,t)) vs time. Note the 
different scales, (e) Second-order pair correlation function (g'^-* (k, — k, vs time. The dotted line shows the idealized 

result (ggL(k,-k,t))^,e = 2 -Hl/(nk(t))^,e for comparison, (f) Number-difference variance (T4(k, — k, vs time, (g) 
Correlation coefficients Ck(i) at |k| = ko vs time. The quantities shown in (d)-(g) are at |k| = ko and the time is in units of 
to = 1/(X\/ Po(0)) = 0.2 ms in all cases. The number of stochastic trajectories used for ensemble averaging is 2000. 
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Thus, in the present example (with the results shown 
in Fig. [SI column 1) the peak density is po{0) = 5 x 10^^ 
m^^, while the TF radius is Rtf — 1-16 x 10"^ m, 
which would correspond to a molecule-molecule scat- 
tering length of floo — 2 nm in a relatively weak trap 
with oscillation frequency uj/2tt = 8 Hz. This gives 
the total initial number of molecules A'o(O) = 1.3 x 10*^. 
The simulations are performed on the spatial domain of 
L = 2.76 X 10"^ m in each dimension, which is twice the 
size of the earlier uniform system L„. The lattice grid 
contained 128^ points. Thus, the lattice spacing in mo- 
mentum space is now Afc = 27r/L ~ 2.28 x 10^ m~^ (half 
that in the uniform system), while the maximum cutoff 
momentum is the same fcmax — 1-46 x 10^ m^^. 

From Fig. O (column 1) we see that the atomic density 
in momentum space (b) and the total number of atoms 
(c) follow closely the predictions of the size-matched uni- 
form system shown in Figs. |2Ib) and ^a). At the 
same time, the mode populations are about eight times 
smaller, which is simply due to the fact that the sim- 
ulation of the nonuniform system is performed within 
a spatial domain which is twice as large, L — 2L,(. 
Accordingly, the respective counting volumes (Afc)^ = 
(27r/L)3 and (Afc„)3 = (27r/i„)3 are different by a 
factor of eight, so that the approximately equal densi- 
ties n(k, t) ~ n„(k, i) give different mode populations 
nk{t) = n(k,t)(Afc)3 andn'^\t) = ri„(k, t)(AA;„)3. Scal- 
ing the uniform result for n„(k, t) with respect to the 
counting volume used in the nonuniform system gives 
n„(k, t)(Afc)^ that approximates well the actual calcu- 
lated value of nk(i) shown in Fig. O 

While the results for the momentum-space density and 
the atom number are in good quantitative agreement 
with the respective results of the size-matched uniform 
system, the same is not true for the correlation func- 
tion g(^)(k, — k, i). The angle-averaged pair correlation 
((/'•^■' (k, — k, at |k| = kg in the present nonuniform 

system is shown in Fig.|5Je) by the full line. Before com- 
paring this result with that of the uniform system [given 

by the idealized result of (7max(k, — k, <) — 2 + l/nk(t)], 
we note that the correct interpretation of the normal- 
ized pair correlation is in the excess of the probability 
of joint detection of pairs of atoms at k and — k rela- 
tive to that probability in an uncorrelated state. Op- 
erationally, this corresponds to determining the number 
of particles in a "detection" volume (Afc)^. Therefore, 
the pair correlation (t/'^^^ (k, — k, i))^^^ for the nonuni- 
form system can be compared with the idealized result 

5max(k, — k, t) = 2 + l/nk(i) provided that nk(t) is eval- 
uated for the same counting volume (Afc)^ as in the 
nonuniform simulation. 

(2) 

The idealized result (7inax(k, — k, t) obtained in this way 
is shown in Fig. (SJe) by the dotted line. We immedi- 
ately see that the strength of the pair correlation for 
the nonuniform system (^(^'(k, — k, t))^_Q is substantially 

lower than (;max(k, — k, i). For comparison, the same 
quantity calculated numerically in the uniform system 



was well approximated by the idealized result [see discus- 
sions of Figs.|3Jb) and^, where it was almost invisible 
behind the solid line for (g(^)(k, — k, 

We note that the correlation between the atoms with 
opposite momenta is still stronger than in an uncorre- 
lated state, ((;(-^^(k, — k, > 1- We have also calcu- 
lated the pair correlation functions ((7^^)(k, k, and 
(g''^-' (k, kgo, As in the uniform case, these were 
given by 2 and 1, respectively, implying that the mode 
mixing has negligible effect on these correlations. 

Another measure of the reduction of the correlation 
strength between the atoms with opposite momenta can 
be obtained using the variance of the number-difference 
fluctuations 14(k, — k, i), Eq. (|31l) . The angle-averaged 
variance {Vyi\<i,—'k.,t))^_g at |k| = fco as a function of 
time t is shown in Fig.[SJf), and is given approximately 
by (l/„(k, -k,i))y,e ^ 0.72, where v = (Afc)^ is the ele- 
mentary volume element of the computational grid. This 
corresponds to 28% squeezing of fluctuations below the 
shot-noise level, which is a significant degradation com- 
pared to the ideal case of 100% squeezing obtained in the 
uniform systems in the absence of molecular depletion or 
s-wave scattering interactions. The reduction of the cor- 
relation strength due to mode mixing is a much stronger 
effect than decorrelation due to the atom-atom recombi- 
nation and s-wave scattering, analyzed in Sees. IIV Bl and 

nvTH 

Finally, in Fig. |3Jg) we plot the correlation coefficient 



Ck(t) 



|(TOk(i))v,e| 



(nk(i))v',e(l + Mt))v,e) 



(36) 



as a function of time t, at |k| — fcg. According to the in- 
equality H24|l . this gives a useful measure of the departure 
[Ck(i) < 1] of the anomalous mode population from the 
maximally correlated, ideal result of Eq. H21() correspond- 
ing to Ck(t) = 1. In addition, the calculated anomalous 
population mk(t) gives us a direct verification of the va- 
lidity of Wick's factorization, according to which the pair 
correlation g^^^ (k, — k, t) in the present example can be 
evaluated as in Eq. (|22|) . We have checked that the off- 
diagonal normal moment (aj^(i)a_k(^)) gives a negligible 
contribution, |(a]j^(i)a_k(0)P/'^k(0 < 1- 



2. Cases 2 and 3: strong inhomogeneity 

Columns 2 and 3 in Fig. show the results of simu- 
lations of a nonuniform system in a tighter trap, with 
w/27r = 19.3 Hz and Rtf = 4.82 x 10"^ m. The total 
initial number of molecules is A'o(O) = 9.4 x 10"^, while all 
other parameters are as in the previous example. Thus, 
in the present case the inhomogeneity in position space 
is stronger. Conversely, the effective coupling g(k) is 
broader in momentum space and therefore the mode mix- 
ing has a stronger effect on the reduction of correlations 
between the atom pairs with opposite momenta. 
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result for g^Ly 
ample, gmLv 



The results presented in column 2 of Fig.[51are obtained 
on the same computational grid as in column 1, i.e., using 
L = 2.76 X 10^^ m and 128^ lattice points. Therefore 
the quantities depending on the counting volume v = 
(Afc)"^ - through the mode occupation number nk(i) = 
n(k, t)(Afc)'^ - are directly comparable. 

Apart from the obvious reduction in the momentum 
space density and total atom number produced as a func- 
tion of time, the fraction of atoms relative to the ini- 
tial number of molecules is almost the same as in the 
example of column 1. The density-density correlation 
function does not reveal significant quantitative change 
either. However, as we explained in the previous subsec- 
tion, the quantitative aspect of the reduction of the corre- 
lation strength should be assessed relative to the idealized 
<(k, — k, t) = 2-1- l/nk(i). In the present ex- 

<(k, — k, t) is off the scale of Fig. EJe), signal- 
ing a much more dramatic reduction in the correlation 
strength than before. This is further seen through the 
number-difference variance, (14(k, — k, i))^^^ — 0.98, im- 
plying only 2% of squeezing below the shot noise-level. 
Similarly, the correlation coefficient Ck(i) is much lower 
than the ideal result of Ck(t) = 1. 

In column 3, Fig.jS] we show the results of simulation 
of the same physical system as in column 2 except that 
the simulation is performed on a smaller (64'^) compu- 
tational grid, with half the length L — 1.38 x 10~^ m. 
Thus, the elementary volume element v = (A/c)^, where 
Afc = 27r/L, is eight times larger and therefore the re- 
sults depending on the counting volume v are scaled by 
a factor of eight. For example, the number-difference 
variance is now given by {VvCk, —'k,t))^^g ~ 0.84 (16% 
squeezing), demonstrating that the larger counting vol- 
ume results in a stronger correlation signal between the 
particle number-difference fluctuations. 
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FIG. 6: (a) Binned distribution of the angle-averaged occu- 
pation numbers (nK(t))v,e a.t t = to; and (b) the respective 
number-difference variance {V—{K.,—K.,t))^^g at |K| — fco as 
a function of time, for the same physical system as in Fig. 
column 1. The graphs can be compared to those of Figs. 
ISJb) and (f), column 1. The {K}-sublattice of the central 
momenta of the bins is setup to correspond to every second 
lattice point of the original computational grid in each dimen- 
sion, so that the bin volume is u = (2Afc)^. The number of 
stochastic trajectories simulated is 4500. 



volume around K. The variance of fluctuations between 
the number of particles in different bins is defined via 



y-(K,K',i) = l 



{fijiit)) + (flK'W) 



(38) 



3. Binning 

The comparison between cases 2 and 3 as a function 
of the counting volume (A/c)^ is similar to the procedure 
of binning except that it is done implicitly - via the vari- 
ation of the length L of the simulated spatial domain in 
each dimension. Since the minimum acceptable length 
must be larger than the characteristic size of the actual 
physical system, the relationship Afc — 2tt/L puts an up- 
per bound on the volume of the elementary bin that can 
be treated in this implicit way. 

To correlate the signals in larger counting volumes, the 
binning must be done explicitly. We introduce the binned 
number operator 



"-k(0 = 




which is defined on a {K}-sublattice corresponding to 
the central momenta of the bins, where iJ(K) is the bin 



where the integration over the bin volume must precede 
the ensemble averaging. 

In Fig. we present the results of simulations of the 
same physical system as in Fig.El column 1, except that 
the results are binned. The bin counting volume for cal- 
culating the occupation numbers n-K{t) — {n-K(t)) and 
the number-difference variance 14r(K,— K, i), Eq. (|38|1 is 
V = (2Afc)'^. This is eight times larger than the elemen- 
tary volume element of the computational grid (Afc)^. 

As expected, the binned variance shows a stronger cor- 
relation signal. The degree of squeezing increases to 
~ 60% ((y-(K,-K,t))^,e ^ 0.4, for |K| = ko), com- 
pared to 28% squeezing in the respective unbinnned vari- 
ance of Fig.[3[f), column 1. Similar results are obtained 
for the cases represented in columns 2 and 3 of Fig. |S1 

In practice, similar issues arise in time-of-flight images 
due to the finite resolution of the imaging system, and 
binning was employed in the recent spatial correlation 
measurements of Ref. . The correlation signal between 
fluctuations of the binned particle numbers was indeed 
found to be stronger for larger bins. 
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V. SUMMARY 

To summarize, we have performed first-principles 
quantum dynamical simulations of dissociation of a BEC 
of molecular dimers into correlated atom pairs in three 
dimensions. The simulations are done using the positive- 
P representation method which allows us to obtain ex- 
act results for the atomic and molecular populations and 
atom-atom correlations in momentum space. We have 
analyzed the effects of molecular depletion, s-wave scat- 
tering interactions, and mode mixing in nonuniform con- 
densates. 

We find that the most useful measure of the strength of 
atom-atom correlations is given by the variance of atom 
number-difference fluctuations in a certain range of mo- 
menta around the central momenta of interest k and k'. 
The strongest correlation signal resulting in squeezing 
of particle number-difference fluctuations below the shot 
noise level is obtained for pairs of equal but opposite mo- 
menta, k' = — k. The degree of squeezing depends on the 
geometry of the system and is stronger for BEC samples 
that are spatially larger and more uniform. 

The major source of atom-atom decorrelation com- 
pared to the ideal result (achievable in uniform systems 
and in the absence of molecular depletion) is mode mixing 
due to inhomogeneity of the molecular BEC. This sug- 
gests that a preferred experimental strategy to obtain a 
strong correlation signal is to start with a large, low den- 
sity sample. The actual correlation strength depends also 
on the counting volume around the central momenta and 
can be increased by binning the signals into larger bins. 

The next important source of decorrelation is the 
atom-atom recombination, which produces increasingly 
large number of molecules in the initially unpopulated 
noncondensate modes as the dissociation proceeds. In 
our example the fraction of these molecules is about 20% 
of the remaining total number once the overall conversion 
reaches 50%. At this stage, the reduction in number- 
difference squeezing due to recombination is about 20%. 

Finally, the simulated s-wave scattering interactions 
are found to have a less severe effect on atom-atom cor- 
relations at least for time durations corresponding to less 
than 10% conversion. For longer time durations, the sim- 
ulations fail due to the limitations of positive-P method, 
unless the atom-atom scattering length is reduced below 
the typical background value. 
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APPENDIX A: DISCRETE AND CONTINUOUS 
FOURIER TRANSFORMS 

To treat the uniform system in a cubic box of side L 
with periodic boundary conditions, we expand the field 
operators (x, t) in terms of the plane- wave momentum 
modes ak(t): 

^'i(x,i) = ;^5]^ak(t)e*-, (Al) 

The inverse transforms are 

= 1571/^'^^ (A3) 
ai{t) = j^ J^d^^{^,t)e^''■-. (A4) 

The mode creation and annihilation operators satisfy the 
usual commutation relation [a-k{t) , al^, {t)] = i5k,k', where 
(5k, k' is the Kroenecker delta function, k = {kx,ky,kz) 
is the momentum (in wave-number units), and ki — 
{2tt /L)ni {i = X, y, z, m 0, ±1, ±2, . . .). 

The treatment of the infinite system in free space 
corresponds to taking the limit of L ^ 00 (A/c = 
2ti/L 0), together with converting the sums = 
En = E„. E„„ into integrals according to X;„. ^ 
^ drii = (L/2iT)Jdki. The corresponding continuous 
Fourier transforms are defined according to 

*i(x,t) = J dkaik,t)e^--, (A5) 

*I (x, t) = J dk at(k, Oe-^*^ ■^ (A6) 

with the inverse transforms of 

a(k, t) = j^l^ J dx *i(x, Oe-^'^ ■^ (A7) 

at(k,t) ^ [a(k,i)]t = J dx*I(x,i)e*-, 

(A8) 

«(-k,t) = J dx^l(x,^)e*''■^ (A9) 

at(-k,t) ^ [a(-k,t)]^ = J dxvl'l(x,i)e-^'^-. 

(AlO) 

The Fourier components a(k) and a^^(k, t) satisfy the 
commutation relation [d(k), a'l'(k')] = (5(k — k'). 

When represented on a discrete computational lattice, 
the continuous Fourier components are approximated by 

a(k,^) = ak(^)/(Afc)^/^ (All) 
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and the continuous delta function by 

^(k-k') = <5k,k'/(Afc)', 



(A12) 



where Afc = 27r/i is the mode spacing in each dimension. 

The Fourier components a(k) (unlike the mode annihila- 
tion operators dk) have units; n{k,t) = (d^(k, <)a(k, t)) 
corresponds to the particle number density in momen- 
tum space, while nk(t) = (dj^(t)dk(t)) is the number of 
particles in the mode k. 



APPENDIX B: FOURIER TRANSFORMS OF 
THE STOCHASTIC FIELDS 

Here, we outline the correspondences between the 
Fourier transforms of the positivc-P stochastic fields and 
those for the field operators. This is important for correct 
implementation of the routines aimed at the calculation 
of various operator moments in Fourier space. 

We recall that the two independent stochastic (c num- 
ber) fields ^'i(x, t) and $i(x, t) are chosen here to corre- 
spond to the annihilation and creation operators 5'i(x, t) 
and ^|(x, according to the following operator corre- 
spondences: 



*i(x,t)^*i(x, t), 
«'I(x,t)^$i(x,t), 



(Bl) 
(B2) 



with similar relationships valid for the molecular field. 
The Fourier transforms of the stochastic fields give: 



a(k,f) = J-x[*i(x,t)](k) 
1 



(27r)3/2 
1 

(27r)3/2 



dx4'i(x,f)e-*-" (B3) 
dx*i(x,i)e-*-^ = a(k,i), 



/3(k,i)=;r^[$i(x,i)](k) 

= J2^2 j d^^^i^,t)e-'''- (B4) 
^ j rfx*I(x,i)e-*- = at(-k,t). 



Thus, the Fourier component /3(k, t) corresponds to 
a^{—\i,t), rather than to d^(k, as might have been ex- 
pected. Therefore, care must be taken in assigning the 
operator correspondences in Fourier space. 

An alternative choice of the positive-P correspon- 
dences between the operators and the stochastic fields 
is 

*i(x,i)^*^(x,t), (B5) 
^l(x,f)^$t(x,t). (B6) 

These stochastic fields are simply the complex conjugates 
of the previous ones: \E'i(x, = '^\{yi,t) and <l>i(x, <) = 
$J(x, t). Their Fourier transforms give: 

[5(k,ir = {^x[$i(x,i)](k)}* 



(B7) 



[/3(k,t)]* = {j-.[$i(x,i)](k)}* 



ik-: 



,zk-x 



(B8) 



it(k,f). 



Thus, the Fourier component d^(k, can be obtained 

through [/3(k, t)]* or /?(— k, t), while d(— k, t) corresponds 
to [5!(k, i)]* and a(— k, f). Using these correspondences, 
one can get access to various operator moments involv- 
ing the Fourier components d(±k, i) and d^(±k, t). For 
example, the normal and anomalous densities can be ob- 
tained via 

(at(k,t)a(k,i)) = ([^(k,ira(k,i)),t 

= (/3(-k,t)a(k,t)),t, (B9) 

(a(k,i)a(-k,i)) = (a(k,t)[5(k,i)]*)st 

= (a(k,i)a(-k,t))3t. (BIO) 
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